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Abstract 

We introduce a family of two-dimensional lattice models of quasicrystals, using a range 
of square hard cores together with a soft interaction based on an aperiodic tiling set. 
Along a low temperature isotherm we find, by Monte Carlo simulation, a first order 
phase transition between disordered and quasicrystalline phases. 
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1. Introduction 

Certain alloys have been shown to exhibit exotic structures in the sense that 
the structures produce sharp X-ray Bragg peaks with symmetries which cannot arise 
from any crystal [1]. We model these materials with equilibrium statistical mechanics 
using a technique due to Levine and Steinhardt [2] which employs an energy ground 
state associated with an "aperiodic tile set" , that is, a set of several geometric shapes 
which only achieve their densest packings in nonperiodic ways, the best known two- 
dimensional example being the kite and dart shapes of Penrose [3]. There is very 
little direct evidence that any equilibrium statistical mechanics model would show a 
well-defined quasicrystalline phase at positive temperature; we know only one model 
for which there is (simulation) evidence [4,5] of such a phase. The theory of solid/fluid 
phase transitions often emphasizes the role of crystal symmetry, going back at least to 
Landau [6,7], and it is of interest to understand how the quasicrystalline state, with 
its exotic form of symmetry [8], fits into this picture. In this regard it is noteworthy 
that the simulations in [4,5] indicate a continuous (high order) transition between the 
fluid and quasicrystalline phases. We revisit this issue here with slightly more realistic 
models and find clear simulation evidence of a discontinuous (first order) transition. 
We contrast our results with related models of hard squares and of Widom-Rowlinson 
type. 

2. The Model 

Our models are variants of that used in [4,5], which we describe first. The lattice 
model of [4,5] uses 16 types of particles, each corresponding to an element of a certain 
subset S C {1,2,3,4,5,6}^ of cardinality 16, associated with an aperiodic tiling set 
due to Ammann [9]. The model consists of configurations of particles on the square 
lattice Z^. It is convenient to interpret each element of iS as a unit-edge square tile, 
with its 4 edges colored using 6 possible colors {1, 2, 3, 4, 5, 6} in one of the 16 groups 
of 4 indicated by S; see Fig. 1. 



1 

2 1 
2 




4 

1 1 
6 




4 

1 2 
5 




6 

1 2 
3 




5 

2 2 
3 




1 

6 4 
1 




2 

5 4 
1 




2 

3 6 
1 


















2 

3 5 
2 




3 

5 3 
5 




3 

5 5 
6 




5 

6 3 
5 




5 

6 5 
6 




5 

6 6 
4 




6 

4 5 
6 




6 

4 6 
4 



Figure 1. Ammann's 16 Wang tiles. 

The model in [4,5] requires precisely one particle/tile centered at each lattice site, 
and each particle interacts with each of its 4 nearest neighbors only, with interaction 
energy or —1, with the negative value when the edge colors match. A canonical 
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ensemble is used, but with density fixed at 1. The model then exhibits a continuous 

phase transition at positive temperature [4,5], with the low temperature phase show- 
ing quasicrystalline structure, including the symmetry [8] it inherits from Ammann's 
aperiodic tiling set. 

We consider variants of the above model, investigating the consequences of two 
changes: extending the size of the square tiles to integer edge length w > 1 (they are 
still "hard" - no overlap of tiles is allowed); and allowing a full range from to 1 of 
volume fraction for the particle configurations, in place of the fixed volume fraction 
of 1 used in [4,5]. Note that the latter adds a second thermodynamic variable to the 
model. 

The description of the interaction is a bit more complicated than in the w = 1 
model of [4,5], since tiles which are close but not touching may still interact, though 
with lower amplitude. To define the interaction we need some notation. A tile-state 
A is an element of 7? x S, with the first coordinate of A representing its center on the 
lattice, and the second coordinate representing its edge colors. We define a pairwise 
interaction H between tiles as follows. For each tile-state A, consider the four "nearest 

neigiiuui i^empLy j aquaiea ui iiif saiiie siZ/t; as J-i. wiiicii suiiuuiiu ^ uixi uu iiUL uvCriap 

it, Isibelled N, S, W and E as in Fig. 2. If another tile-state A' is disjoint with A but 





N 












w 


A * 


E 




S 





Figure 2. The interaction energy depends on the 
£^ distance between particle centers (dashed line). 

overlaps one of these four squares - it cannot possibly overlap more than one - its 
energy of interaction with A, H{A, A'), is as follows. Assume without loss of generality 
that A' overlaps N. Then H{A, A'): is negative if and only if the north edge of A has 
the same color as the south edge of A', and positive otherwise; and has absolute value 
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given by \H(A, A')\ = {3w — d)/{2w), where d is the distance between the centers 
of A and A' . (The £^ distance between (x, and {x',y') is \x — x'\ + \y — y'\-) If A' 
overlaps 5", W, or E, then the energy H is defined similarly. 

Given an inverse temperature f3, a chemical potential ^ (common to all tiles/particle 
types), and a finite volume V C Z^, we consider the grand canonical ensemble: the 
probability Vp^^ on configurations C of tile-states in V which is given by 

e^p[-miC)-fiAr{C))] 

yi3,AC) = , (Ij 

where M{C) is the number of tiles-states in C, 'H(C) = J^a^A'eC -^(^' total 
energy of C, and Z = ^(/U, /3) is the appropriate normalization. In our simulations we 




Figure 3: A boundary shell (with edges colored instead of numbered). 

use "square" volumes V = fl [— /c, /c]^ with k an integer multiple of w. We use 
boundary conditions as follows. Starting from a perfect tiling T by tile-states, that 
is, a configuration of tiles covering in which all edge colors match, we remove all 
the tiles except those just outside V; this leaves a (fixed) boundary "shell" ; see Fig. 3. 
(We use 10,000 different boundary shells arising from random translations of a single 
tiling.) 

It remains to define the order parameter, x- Each tile-state has a "long/short 
pattern" given by its 4 edge colors, where we call colors 1 and 2 "short" and colors 3, 
4, 5 and 6 "long". Note that, for any tile, horizontal edges are either both short or 
both long, and vertical edges are either both short or both long, so there are 4 possible 
long/short patterns for a tile. Any tile-state A centered "near" the site {wi,wj) G V, 
namely centered in [wi — w/2,wi + w/2) x [wj — w/2,wj + w/2), will be said to agree 
with T if it has the same long/short pattern as the tile-state in T centered at {wi, wj). 
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Now, consider the thick "shells" Vj = V — Given a configuration C in V 

and boundary from T, let j* be the largest integer such that at least 80% of sites 
{wi, wj) in Vj* have a tile-state centered nearby (in the above sense) which agrees 
with T. The order parameter x is defined as the corresponding normalized volume, 

X{C) = \V,*\/\V\. 

3. Simulations, and Results 

Wc ran Markov-chain Monte Carlo simulations of the model with tile widths w = 
1, 2, and 3. For each tile width w we simulated volumes |V| = (20to)^, (40w)^, (60u;)^, 
{80w)^ and (lOOw)^. (We will also use the notation V = \V\/w'^, so V represents the 
maximum possible number of tile-states in a single configuration.) The Monte Carlo 
steps were designed to produce the target distribution of equation (1). The basic Monte 
Carlo step involves a choice between four types of moves: moving a particle (locally), 
changing a particle type, removing a particle, and adding a particle. As usual in 
grand canonical Monte Carlo simulation, given a configuration C a trial configuration 
C is introduced, effecting one of the moves described above, and is accepted with 
probability determined by the desired limiting distribution, that is, equation (1). 

More precisely, the basic Monte Carlo step is as follows. Let C{n) be the con- 
figuration at the nth step of the Monte Carlo chain, and choose a random lattice 
site X. Select k G {1,2,3,4} with probability p^, where the pk are the probabilities 
of attempting the four types of Monte Carlo moves, with the pk summing to 1 and 
P3 = Pa- (We choose the pk somewhat arbitrarily.) Without loss of generality suppose 
X & R := [wi — w/2, wi + w/2) x [wj — w /2,wj + w /2). First suppose k G {1, 2}. If 
there is no tile-state centered in R, then C{n -\- 1) = C{n); otherwise, let A be the 
(unique) tile-state centered in R, and produce a trial configuration C{n)' by: if /c = 1, 
moving the center of to a site in R chosen uniformly at random; if /c = 2, replacing 
A with a different tile-state centered at the same lattice site. Now suppose A; = 3. If 
there is a tile-state A centered in i?, then C{n)' is obtained by removing A; otherwise 
C{n + 1) = Cin). Finally suppose k = A. If there is no tile-state centered in i?, 
then let Cin)' be the configuration obtained from C(n) by adding a tile-state chosen 
uniformly at random and centered at x; otherwise C(n -|- 1) = C{n). In cases where 
a trial configuration C{n)' is introduced, we take C{n+ 1) = C{n)' with probability 
Q = min(l, q); otherwise we take C{n+1) = C{n). Here q is given by 

q := min{l,e'5t^(^)-^(^')+'^*(-^(^')-^(^))]}, (2) 

where we have written C for C{n) and C for C(n)', and take 'H{C') = oo for trial 
configurations C with overlapping tiles. Here fi* = + /3~^ ln(16s^) is chosen so that 
the Monte Carlo configurations C(n) have limiting probability distribution given by 
equation (1). (The value //* arises because of the following. Starting with a configura- 
tion C, the probability of proposing a trial configuration C which removes tile-state 
A from C is 16s^ times the probability that, starting with the configuration C, a trial 
configuration C" = C is proposed which adds the tile-state A to C .) 

For w — 2 and w = 3 we find clear evidence of discontinuous (first order) phase 
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Figure 4: Data for w = 2 and /3 = 1.5. 

transitions. At the transitions (at /3 = 1.5, near /j, = 1.4 for w = 2 and near fi = 2.0 for 
w = 3) the average volume fraction (p and energy per volume Ti/V exhibit developing 
jump discontinuities as V increases; see Figs. 4-5. Furthermore the (volumetric) heat 
capacity, Cy (defined as the derivative of Ti/V with respect to T = 1//3), exhibits a 
developing delta function as V increases (see Figs. 4-5). 
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Figure 5: Data for w = 3 and (3 = 1.5. 



From our data in Fig. 6, for the w = 1 model, there appears to be a continuous 
transition (at /3 — 1.5 near fj, = 0.8), from a disordered fluid phase to a phase with 
quasicrystalline order, similar to the transition in [4,5] of the canonical ensemble at 
density 1. 
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Figure 6: Data for w = 1 and /3 = 1.5. 

To quantify the developing delta function in Cy noted above for = 2, 3 we 
measure, from Figs. 4 and 5, the maximum value h of Cy divided by the "width" of the 
graph of Cv at h/2 (see Fig. 7). The order parameter x also develops a discontinuity 
at the transition, signaling the onset of quasicrystalline symmetry: see Figs. 4-5. 
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Figure 7. Evidence of Cy developing into a delta function for a) w = 2, and b) w = 3. 
4. Comparisons 

We next compare our results with related models, in particular hard squares and 
Widom-Rowlinson models on 1?; see [10] for a review of, and earlier references to, the 
former and [11] for a review of, and earlier references to, the latter. 

At sufficiently high temperature the soft interaction between our particles must 
become negligible compared to the hard core so we first review what is known about 
hard square systems. 

For hard squares on 1? the energy is zero for all configurations so it is convenient 
to use a grand canonical ensemble with fixed temperature 1 and variable chemical 
potential /i. Hard squares with edge size w — 1 ("point hard core") do not have a phase 
transition as /i is varied, since the variables at each site, 'occupied' or 'unoccupied', 
are independent. For fixed w > 2 the densest configurations, of volume fraction 1 and 
corresponding to /U — ?■ oo, are degenerate, allowing parallel rows or columns of squares 
to slide independently in what is called a "columnar" phase. It is generally accepted 
from simulation [10] that the model with w = 2 has a continuous transition at volume 
fraction about 0.92. We only know of one paper on the model with w — 3 [10]; the 
authors claim a first order freezing transition beginning at a density above 0.91, but 
are unable to show developing discontinuities in any thermodynamic quantities. 

We now compare our model to hard squares on Z^. Consider a version of our 
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quasicrystal model with w >2 but in a canonical ensemble with volume fraction fixed 
at 1. Based on the results of [4,5], we expect this model to have a columnar phase at 
high temperature, and a quasicrystalline phase at low temperature, with a transition 
between them at some unknown temperature Tg. Now considering the same model 
except at lower volume fraction or chemical potential, one expects the transition to the 
quasicrystalline phase to move to lower temperature since the soft interaction should 
have less effect as the density is lowered. We conclude that above temperature Tc our 
model should behave qualitatively like hard squares, with disordered and columnar 
phases, but without a phase with quasicrystalline order. 

Absent from the above discussion is the effect on entropy, in our model, of the 
multiple types of particle. This effect is emphasized in Widom-Rowlinson models, 
which we now discuss. The multitype Widom-Rowlinson model with q colors, on 7?, 
is defined by the hard core condition that nearest neighbor sites cannot be occupied by 
different colors. The main conclusion of such work is the existence (for large enough q) 
of two types of high density phase transitions: one transition at which spatial symmetry 
is broken, producing different population densities for the two sublattices, but without 
breaking the symmetry between particle types, and a second (discontinuous) transition 
at higher density at which the symmetry between particle types is broken, so that one 
particle type dominates the density [11]. These models, with diamond-shaped instead 
of square hard cores, are not simply related to our models, but there are two interesting 
variants treated by Georgii and Zagrebnov [12] which are closely related. The first 
model has a square-shaped hard core, with square width w = 2, between particles 
of different color; the second model adds to this a smaller diamond-shaped hard core 
between nearest neighbor particles of the same color. One could think of these models 
as employing a hard core of size w = 2 between all particles, with a "cancelling" (or 
attractive) interaction between particles of like color on sites at separation \/2 for the 
second model, and for sites at separation 1 and for the first model. Thus the models 
can be thought of as hard squares with an added particle-type-dependent short range 
attraction, somewhat like our quasicrystal models. Both Widom-Rowlinson models 
show only one transition, which is discontinuous. For the first model the high density 
phase breaks particle type symmetry but not spatial symmetry, and for the second 
model the high density phase breaks the symmetry of both space and particle type. 
In particular the second model seems close to our {w — 2) quasicrystal model, and our 
results bear this out, with a discontinuous transition corresponding to broken spatial 
and particle-type symmetries. The correspondence might be closer if the attraction 
and repulsion in our model was a hard rather than soft interaction, as in [13]. 

It would have been useful to simulate our model over a range of temperatures, 
but we found this prohibitively expensive in computation time. For this reason we 
investigated a technique noted in [14], in which the initial state for the simulation is 
an energy ground state. However this technique proved to be unreliable here unless 
the runs were of comparable length to those starting from vacuum, thus yielding no 
advantage (see Fig. 8). (We suspect that simulations starting from an energy ground 
state get restricted to a narrow portion of phase space more easily than those starting 
from the vacuum.) 
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Figure 8. Comparison of order parameter values of long and short runs starting from 
both vacuum and perfect tiling, for w — 1 and V — 100^. (The long runs are 100 times 
longer than the short runs.) 

5. Conclusion 

We have introduced an extension of a two dimensional lattice gas model of a 
quasicrystal [4,5]. That model was studied as a function of temperature at fixed density 
1, with simulations showing a continuous phase transition between a disordered (fluid) 
phase and a phase with quasicrystalline order. When we introduce vacancies into the 
model we find that the basic character does not seem to change: there still appears to 
be only the two phases, disordered and quasicrystalline, with a continuous transition 
between them. 

We then generalized the model by extending the range of the hard core as well as 
allowing variation in the density, and found, by simulation along a (low temperature) 
isotherm, that the model then exhibits a discontinuous transition between a low density 
disordered and high density quasicrystalline phase. The high density quasicrystalline 
phase breaks symmetries of space (a square sublattice is preferentially populated) and 
particle type (the particles inherit the quasicrystalline boundary structure). Along an 
isotherm at sufficiently high temperature we expect the model to exhibit the transition 
of simple hard squares (continuous for tx; = 2 and possibly discontinuous for w = 3), 
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between the low density disordered and high density columnar phases. At fixed high 
density or chemical potential we expect a transition from the quasicrystalline phase to 
a columnar phase as temperature is raised from zero. 

This family of models with unusual symmetry adds a special flavor to the more 
traditional lattice gas models with hard core. Its complicated multiparticle structure 
seems at present to defy proof of a transition, and is also expensive to simulate. Further 
progress in filling out the phase structure of the model would be highly desirable. 
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